home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Over 1,000 Windows 95 Programs
/
Over 1000 Windows 95 Programs (Microforum) (Disc 1).iso
/
1249
/
root.t
< prev
next >
Wrap
Text File
|
1997-04-18
|
34KB
|
1,681 lines
%
% "root.t" demonstates techniques for the numerical
% solution of non-linear equations
%
% Sample program for the T Interpreter by:
%
% Stephen R. Schmitt
% 962 Depot Road
% Boxborough, MA 01719
%
const MATVECT_EPSILON : real := 1.0E-15
const TOLERANCE : real := 0.00000001
const MAX_ITER : int := 50
const DIM : int := 10
type rvector : array[DIM] of real
type ivector : array[DIM] of int
type rmatrix : array[DIM,DIM] of real
type complex : record
re : real
im : real
end record
type cvector : array[DIM] of complex
program
do_bisection
do_newton_approx
do_newton_exact
do_richmond_approx
do_richmond_exact
do_combined
do_brent
do_deflation
do_lin_bairstow
do_newton_two_fncs
do_newton_multi_roots
do_newton_snle
end program
%---------------------------------------------------------------------
% demonstrations of the numerical solution of non-linear equations
%---------------------------------------------------------------------
%
% Fx returns the value of: f(x) = exp(x) - 3*x^2
%
% d_Fx returns the first derivative: f'(x) = exp(x) - 6*x
%
% dd_Fx returns the second derivative: f"(x) = exp(x) - 6
%
function Fx( x : real ) : real
return exp(x) - 3*x^2
end function
function d_Fx( x : real ) : real
return exp(x) - 6*x
end function
function dd_Fx( x : real ) : real
return exp(x) - 6
end function
const Fx_msg : string := "solving: f(x) = exp(x) - 3*x^2"
%
% returns value of either of two sumultaneous non-linear equations:
%
% f0(x0,x1) = x1*exp(x0) - 2
% f1(x0,x1) = x0^2 + x1 - 4
%
function Snle( var x : rvector, index : int ) : real
case index of
value 0:
return x[1] * exp( x[0] ) - 2
value 1:
return x[0]^2 + x[1] - 4
value:
return 1
end case
end function
const f0_msg : string := "f0(x0,x1) = x1*exp(x0) - 2"
const f1_msg : string := "f1(x0,x1) = x0^2 + x1 - 4"
%
% F1xy returns the value of: f1(x,y) = x^2 + y^2 - 1
%
% F2xy returns the value of: f2(x,y) = x^2 - y^2 + 0.5
%
function F1xy( x, y : real ) : real
return x^2 + y^2 - 1
end function
function F2xy( x, y : real ) : real
return x^2 - y^2 + 0.5
end function
const F1xy_msg : string := "f1(x,y) = x^2 + y^2 - 1"
const F2xy_msg : string := "f2(x,y) = x^2 - y^2 + 0.5"
%
% test Bisection
%
procedure do_bisection
var low, high, root : real
low := 2
high := 4
put "bisection method\n"
put Fx_msg
put "guess interval is: ", low, "...", high
if Bisection( low, high, root ) then
put "root = ", root
else
put "could not solve"
end if
end procedure
%
% test Newton_approx
%
procedure do_newton_approx
var root : real := 3
var num_roots, i : int
put "\nNewton's method for an approximate derivative\n"
put Fx_msg
put "initial guess is: ", root
if Newton_approx( root ) then
put "root = ", root
else
put "could not solve"
end if
end procedure
%
% test Newton_exact
%
procedure do_newton_exact
var root : real := 3
put "\nNewton's method for a user provided derivative\n"
put Fx_msg
put "initial guess is: ", root
if Newton_exact( root ) then
put "root = ", root
else
put "could not solve"
end if
end procedure
%
% test Richmond_approx
%
procedure do_richmond_approx
var root : real := 3
put "\nRichmond's method with an approximate derivative\n"
put Fx_msg
put "initial guess is: ", root
if Richmond_approx( root ) then
put "root = ", root
else
put "could not solve"
end if
end procedure
%
% test Richmond_exact
%
procedure do_richmond_exact
var root : real := 3
put "\nRichmond's method with a user provided derivative\n"
put Fx_msg
put "initial guess is: ", root
if Richmond_exact( root ) then
put "root = ", root
else
put "could not solve"
end if
end procedure
%
% test Combined
%
procedure do_combined
var low, high, root : real
low := 2
high := 4
put "\nCombined method\n"
put Fx_msg
put "guess interval is: ", low, "...", high
if Combined( low, high, root ) then
put "root = ", root
else
put "could not solve"
end if
end procedure
%
% test Brent
%
procedure do_brent
var low, high, root : real
low := 2
high := 4
put "\nBrent's method\n"
put Fx_msg
put "guess interval is: ", low, "...", high
if Brent( low, high, root ) then
put "root = ", root
else
put "could not solve"
end if
end procedure
%
% test Deflate_polynomial
%
procedure do_deflation
const SMALL : real := 0.00000001
var coeff, roots : rvector
var poly_order, num_roots, i : int
poly_order := 3
coeff[0] := -40
coeff[1] := 62
coeff[2] := -23
coeff[3] := 1
put "\nDeflating Polynomial method\n"
put "the polynomial:"
for decreasing i := poly_order...0 do
if i > 1 then
put coeff[i], "*X^", i, " + "...
elsif i = 1 then
put coeff[i], "*X + "...
elsif i = 0 then
put coeff[i]...
end if
end for
put " = 0"
if Deflate_polynomial( coeff, 1.1, roots,
num_roots, poly_order ) then
put "has roots:"
for i := 0...poly_order-1 do
put "root ", i + 1, " = ", roots[i]
end for
else
put "could not solve"
end if
end procedure
%
% test Lin_Bairstow
%
procedure do_lin_bairstow
var coeff : rvector
var roots : cvector
var r : complex
var num_roots, poly_order, i : int
poly_order := 3
coeff[0] := 8
coeff[1] := 4
coeff[2] := 2
coeff[3] := 1
put "\nLin-Bairstow method\n"
put "the polynomial:"
for decreasing i := poly_order...0 do
if i > 1 then
put coeff[i], "*X^", i, " + "...
elsif i = 1 then
put coeff[i], "*X + "...
elsif i = 0 then
put coeff[i]...
end if
end for
put " = 0"
if Lin_Bairstow( coeff, roots, poly_order ) then
put "has roots:"
for i := 0...poly_order-1 do
r.re := roots[i].re
r.im := roots[i].im
put cstr( r )
end for
else
put "could not solve"
end if
end procedure
%
% test Newton2functions
%
procedure do_newton_two_fncs
var rootX, rootY : real
rootX := 1
rootY := 3
put "\nNewton's method for two equations\n"
put "solving: "
put F1xy_msg
put F2xy_msg
put "initial guess for x is ", rootX
put "initial guess for y is ", rootY
if Newton2functions( rootX, rootY ) then
put "root X = ", rootX
put "root Y = ", rootY
else
put "could not solve"
end if
end procedure
%
% test NewtonMultiRoots
%
procedure do_newton_multi_roots
var roots : rvector
var num_roots, i : int
roots[0] := 2
put "\nNewton's method for multiple roots\n"
put Fx_msg
put "initial guess is: ", roots[0]
if NewtonMultiRoots( roots, num_roots, DIM ) then
for i := 0...num_roots-1 do
put "root ", i + 1, " = ", roots[i]
end for
else
put "could not solve"
end if
end procedure
%
% test NewtonSimNLE
%
procedure do_newton_snle
var Xguess : rvector
var num_eqns, i : int
num_eqns := 2
Xguess[0] := -0.16
Xguess[1] := 2.7
put "\nNewton's method for simultaneous non-linear equations\n"
put "solving for:"
put f0_msg
put f1_msg
for i := 0...num_eqns-1 do
put "initial guess for X", i, " = ", Xguess[i]
end for
if NewtonSimNLE( Xguess, num_eqns ) then
put "roots are:"
for i := 0...num_eqns-1 do
put "X", i, " = ", Xguess[i]
end for
else
put "could not solve"
end if
end procedure
%---------------------------------------------------------------------
% functions for numerical solution of non-linear equations
%---------------------------------------------------------------------
%
% Solve for a root of a function of a single variable by selecting
% an initial interval that contains the root and searching
%
function Bisection( low, high : real, var root : real ) : boolean
if Fx( low ) * Fx( high ) > 0 then
return false
end if
loop
exit when abs( high - low ) < TOLERANCE
% update guess
root := ( low + high ) / 2
if Fx( root ) * Fx( high ) > 0 then
high := root
else
low := root
end if
end loop
return true
end function
%
% Solve for a root of a function of a single variable by selecting
% an initial guess and using the function value and an approximate
% derivative to search for the solution
%
function Newton_approx( var root : real ) : boolean
var iter : int := 0
var h, diff : real
loop
h := 0.01 * root
if abs( root ) < 1 then
h := 0.01
end if
% calculate guess refinement
diff := 2 * h * Fx( root )
diff := diff / ( Fx( root+h ) - Fx( root-h ) )
% update guess
root := root - diff
iter := iter + 1
exit when iter > MAX_ITER or abs( diff ) < TOLERANCE
end loop
if abs( diff ) <= TOLERANCE then
return true
else
return false
end if
end function
%
% Solve for a root of a function of a single variable by selecting
% an initial guess and using the function value and the
% derivative to search for the solution
%
function Newton_exact( var root : real ) : boolean
var iter : int := 0
var diff : real
loop
% calculate guess refinement
diff := Fx( root ) / d_Fx( root )
% update guess
root := root - diff
iter := iter + 1
exit when iter > MAX_ITER or abs( diff ) < TOLERANCE
end loop
if abs( diff ) <= TOLERANCE then
return true
else
return false
end if
end function
%
% Solve for a root of a function of a single variable by selecting
% an initial guess and using the function value and an approximation
% of the first and second derivative to search for the solution
%
function Richmond_approx( var root : real ) : boolean
var iter : int
var h, diff : real
var f1, f2, f3 : real
var fd1, fd2 : real
loop
h := 0.01 * root
if abs( root ) < 1 then
h := 0.01
end if
f1 := Fx( root - h ) % f(x-h)
f2 := Fx( root ) % f(x)
f3 := Fx( root + h ) % f(x+h)
fd1 := ( f3 - f1 ) / ( 2 * h ) % f'(x)
fd2 := ( f3 - 2 * f2 + f1 ) / sqrt( h ) % f"(x)
% calculate guess refinement
diff := f1 * fd1 / ( fd1 ^ 2 - 0.5 * f1 * fd2 )
% update guess
root := root - diff
iter := iter + 1
exit when iter > MAX_ITER or abs( diff ) < TOLERANCE
end loop
if abs( diff ) <= TOLERANCE then
return true
else
return false
end if
end function
%
% Solve for a root of a function of a single variable by selecting
% an initial guess and using the function value and an approximation
% of the first and second derivative to search for the solution
%
function Richmond_exact( var root : real ) : boolean
var iter : int
var diff, f1 : real
var fd1, fd2 : real
iter := 0
loop
f1 := Fx( root )
fd1 := d_Fx( root )
fd2 := dd_Fx( root )
% calculate guess refinement
diff := f1 * fd1 / ( fd1^2 - 0.5 * f1 * fd2 )
% update guess
root := root - diff
iter := iter + 1
exit when iter > MAX_ITER or abs( diff ) < TOLERANCE
end loop
if abs( diff ) <= TOLERANCE then
return true
else
return false
end if
end function
%
% Brent's method uses interpolation in a variation of
% the bi-section method.
%
function Brent( low, high : real, var root : real ) : boolean
const SMALL : real := 0.0000001 % epsilon
const VERY_SMALL : real := 0.0000000001 % near zero
var iter : int
var a, b, c : real
var fa, fb, fc : real
var d, e, tol : real
var small1, small2, small3 : real
var p, q, r : real
var s, xm : real
iter := 1
a := low
b := high
c := high
fa := Fx( low )
fb := Fx( high )
fc := fb
% check that the guesses contain the root
if fa * fb > 0 then
return false
end if
% start loop to refine the guess for the root
loop
exit when iter > MAX_ITER
iter := iter + 1
if fb * fc > 0 then
c := a
fc := fa
e := b - a
d := e
end if
if abs(fc) < abs(fb) then
a := b
b := c
c := a
fa := fb
fb := fc
fc := fa
end if
tol := 2 * SMALL * abs(b) + TOLERANCE / 2
xm := ( c - b ) / 2
if abs( xm ) <= tol or abs( fb ) <= VERY_SMALL then
root := b
return true
end if
if abs( e ) >= tol and abs( fa ) > abs( fb ) then
% perform the inverse quadratic interpolation
s := fb / fa
if abs(a - c) <= VERY_SMALL then
p := 2 * xm * s
q := 1 - s
else
q := fa / fc
r := fb / fc
q := s * (2 * xm * q * (q - r) - (b - a) * (r - 1))
q := (q - 1) * (r - 1) * (s - 1)
end if
% determine if improved guess is inside the range
if p > 0 then
q := -q
end if
p := abs( p )
small1 := 3 * xm * q * abs(tol * q)
small2 := abs(e * q)
if small1 < small2 then
small3 := small1
else
small3 := small2
end if
if 2 * p < small3 then
% interpolation successfull
e := d
d := p / q
else
% use bisection because the
% interpolation did not succeed
d := xm
e := d
end if
else
% use bisection because the range
% is slowly decreasing
d := xm
e := d
end if
% copy most recent guess to variable a
a := b
fa := fb
% evaluate improved guess for the root
if abs(d) > tol then
b := b + d
else
if xm > 0 then
b := b + abs( tol )
else
b := b - abs( tol )
end if
end if
fb := Fx( b )
end loop
return false
end function
%
% This process combines the bi-section and newton techniques
%
function Combined( low, high : real, var root : real ) : boolean
var iter : int
var h, diff : real
iter := 0
loop
iter := iter + 1
h := 0.01 * root
if abs( root ) < 1 then
h := 0.01
end if
% calculate guess refinement
diff := 2 * h * Fx( root ) / ( Fx( root + h ) - Fx( root - h ) )
root := root - diff
% check if Newton's method yields a refined guess
% outside the range (low, high)
if root < low or root > high then
% apply Bisection method for this iteration
root := ( low + high ) / 2
if Fx( root ) * Fx( high ) > 0 then
high := root
else
low := root
end if
end if
exit when iter > MAX_ITER or abs( diff ) < TOLERANCE
end loop
if abs( diff ) <= TOLERANCE then
return true
else
return false
end if
end function
%
% The deflating polynomial method combines a polynomial reduction
% step with Newton's method. As each (real) root is found, the
% order of the polynomial to be solved is reduced.
%
function Deflate_polynomial( var coeff : rvector,
initGuess : real,
var roots : rvector,
numRoots, polyOrder : int ) : boolean
var a b, c : rvector
var diff, z, X : real
var i, iter, n : int
iter := 1
n := polyOrder + 1
X := initGuess
numRoots := 0
for i := 0...n-1 do
a[i] := coeff[i]
end for
loop
exit when iter > MAX_ITER or n <= 1
iter := iter + 1
z := X
b[n-1] := a[n-1]
c[n-1] := a[n-1]
for decreasing i := n-2...1 do
b[i] := a[i] + z * b[i+1]
c[i] := b[i] + z * c[i+1]
end for
b[0] := a[0] + z * b[1]
diff := b[0] / c[1]
X := X - diff
if abs(diff) <= TOLERANCE then
iter := 1 % reset iteration counter
n := n - 1
roots[numRoots] := X
numRoots := numRoots + 1
% update deflated roots
for i := 0...n-1 do
a[i] := b[i + 1]
% get the last root
if n = 2 then
n := n - 1
roots[numRoots] := -a[0]
numRoots := numRoots + 1
end if
end for
end if
end loop
return true
end function
%
% The Lin-Bairstow method improves on the deflation method to solve
% for the complex roots of the following polynomial:
%
% y = coeff(0) + coeff(1) X + coeff(2) X^2 +...+ coeff(n) X^n
%
% parameters:
%
% coeff must be an array with at least order+1 elements.
% roots output array of roots
% order order of polynomial
%
function Lin_Bairstow( var coeff : rvector,
var roots : cvector,
order : int ) : boolean
const SMALL : real := 0.00000001
var a, b, c, d : rvector
var alfa1, alfa2 : real
var beta1, beta2 : real
var delta1, delta2, delta3 : real
var i, j, k : int
var count : int
var n : int
var n1 : int
var n2 : int
n := order
n1 := n + 1
n2 := n + 2
% is the coefficient of the highest term zero?
if abs( coeff[0] ) < SMALL then
return false
end if
for i := 0...n do
a[n1-i] := coeff[i]
end for
% is highest coeff not close to 1?
if abs( a[1]-1 ) > SMALL then
% adjust coefficients because a(1) != 1
for i := 2...n1 do
a[i] := a[i] / a[1]
end for
a[1] := 1
end if
% initialize root counter
count := 0
%
% start the main Lin-Bairstow iteration loop
% initialize the counter and guesses for the
% coefficients of quadratic factor:
%
% p(x) = x^2 + alfa1 * x + beta1
%
loop
alfa1 := 1
beta1 := 1
loop
b[0] := 0
d[0] := 0
b[1] := 1
d[1] := 1
j := 1
k := 0
for i := 2...n1 do
b[i] := a[i] - alfa1 * b[j] - beta1 * b[k]
d[i] := b[i] - alfa1 * d[j] - beta1 * d[k]
j := j + 1
k := k + 1
end for
j := n - 1
k := n - 2
delta1 := d[j]^2 - ( d[n] - b[n] ) * d[k]
alfa2 := ( b[n] * d[j] - b[n1] * d[k] ) / delta1
beta2 := ( b[n1] * d[j] - ( d[n] - b[n] ) * b[n] ) / delta1
alfa1 := alfa1 + alfa2
beta1 := beta1 + beta2
exit when abs( alfa2 ) < TOLERANCE and abs( beta2 ) < TOLERANCE
end loop
delta1 := alfa1^2 - 4 * beta1
if delta1 < 0 then % imaginary roots
delta2 := sqrt( abs( delta1 ) ) / 2
delta3 := -alfa1 / 2
roots[count].re := delta3
roots[count].im := +delta2
roots[count+1].re := delta3
roots[count+1].im := -delta2
else % roots are real
delta2 := sqrt( delta1 )
roots[count].re := ( delta2 - alfa1 ) / 2
roots[count].im := 0
roots[count+1].re := ( delta2 + alfa1 ) / ( -2 )
roots[count+1].im := 0
end if
% update root counter
count := count + 2
% reduce polynomial order
n := n - 2
n1 := n1 - 2
n2 := n2 - 2
% for n >= 2 calculate coefficients of
% the new polynomial
if n >= 2 then
for i := 1...n1 do
a[i] := b[i]
end for
end if
exit when n < 2
end loop
if n = 1 then
% obtain last single real root
roots[count].re := -b[2]
roots[count].im := 0
end if
return true
end function
%
% This function uses Newton's method adapted for the solution of
% two simultaneous non-linear equations.
%
function Newton2functions( var rootX, rootY : real ) : boolean
var Jacob : real
var fx0, fy0 : real
var hX, hY : real
var diffX, diffY : real
var fxy, fxx : real
var fyy, fyx : real
var X : real
var Y : real
var iter : int
X := rootX
Y := rootY
iter := 1
loop
hX := 0.01 * rootX
if abs( rootX ) < 1 then
hX := 0.01
end if
hY := 0.01 * rootY
if abs(rootY) < 1 then
hY := 0.01
end if
fx0 := F1xy( X, Y )
fy0 := F2xy( X, Y )
fxx := ( F1xy( X + hX, Y ) - F1xy( X - hX, Y ) ) / 2 / hX
fyx := ( F2xy( X + hX, Y ) - F2xy( X - hX, Y ) ) / 2 / hX
fxy := ( F1xy( X, Y + hY ) - F1xy( X, Y - hY ) ) / 2 / hY
fyy := ( F2xy( X, Y + hY ) - F2xy( X, Y - hY ) ) / 2 / hY
Jacob := fxx * fyy - fxy * fyx
diffX := ( fx0 * fyy - fy0 * fxy ) / Jacob
diffY := ( fy0 * fxx - fx0 * fyx ) / Jacob
X := X - diffX
Y := Y - diffY
iter := iter + 1
exit when iter > MAX_ITER
exit when abs( diffX ) < TOLERANCE and abs( diffY ) < TOLERANCE
end loop
rootX := X
rootY := Y
if abs( diffX ) <= TOLERANCE and abs( diffY ) <= TOLERANCE then
return true
else
return false
end if
end function
%
% This function uses Newton's method adapted for the solution of
% multiple simultaneous non-linear equations.
%
function NewtonMultiRoots ( var roots : rvector,
var numRoots : int,
maxRoots : int ) : boolean
var iter, i : int
var h, diff : real
var f1, f2, f3 : real
var root : real
numRoots := 0
root := roots[0]
loop
iter := 0
loop
h := 0.01 * root
if abs( root ) < 1 then
h := 0.01
end if
f1 := Fx( root - h )
f2 := Fx( root )
f3 := Fx( root + h )
if numRoots > 0 then
for i := 0...numRoots-1 do
f1 := f1 / ( root - h - roots[i] )
f2 := f2 / ( root - roots[i] )
f3 := f3 / ( root + h - roots[i] )
end for
end if
% calculate guess refinement
diff := 2 * h * f2 / ( f3 - f1 )
% update guess
root := root - diff
iter := iter + 1
exit when iter > MAX_ITER or abs( diff ) <= TOLERANCE
end loop
if abs(diff) <= TOLERANCE then
roots[numRoots] := root
numRoots := numRoots + 1
if root < 0 then
root := 0.95 * root
elsif root > 0 then
root := 1.05 * root
else
root := 0.05
end if
end if
exit when iter > MAX_ITER or numRoots >= maxRoots
end loop
if numRoots > 0 then
return true
else
return false
end if
end function
%
% This function uses another adaptation of Newton's method for the
% solution of multiple simultaneous non-linear equations. It uses
% matrix arithmetic functions LUDecomp and LUBackSubst.
%
function NewtonSimNLE( var X : rvector,
numEqns : int ) : boolean
var Xdash : rvector
var Fvector : rvector
var index : ivector
var Jmat : rmatrix
var i, j : int
var moreIter : boolean
var iter : int
var rowSwapFlag : int
var h : real
var dummy : boolean
iter := 0
loop
iter := iter + 1
% copy the values of array X into array Xdash
for i := 0...numEqns-1 do
Xdash[i] := X[i]
end for
% calculate the array of function values
for i := 0...numEqns-1 do
Fvector[i] := Snle( X, i )
end for
% calculate the Jmat matrix
for i := 0...numEqns-1 do
for j := 0...numEqns-1 do
% calculate increment in variable number j
if abs( X[j] ) > 1 then
h := 0.01 * X[j]
else
h := 0.01
end if
Xdash[j] := Xdash[j] + h
Jmat[i,j] := ( Snle( Xdash, i ) - Fvector[i] ) / h
% restore incremented value
Xdash[j] := X[j]
end for % j
end for % i
% solve for the guess refinement vector
dummy := LUDecomp( Jmat, index, numEqns, rowSwapFlag )
LUBackSubst( Jmat, index, numEqns, Fvector )
% clear the more-iteration flag
moreIter := false
% update guess and test for convergence
for i := 0...numEqns-1 do
X[i] := X[i] - Fvector[i]
if abs( Fvector[i] ) > TOLERANCE then
moreIter := true
end if
end for
% check iteration limit
if moreIter then
if iter > MAX_ITER then
moreIter := false
else
moreIter := true
end if
end if
exit when not moreIter
end loop
return not moreIter
end function
%
% This function performs LU decomposition of a matrix of real values.
%
function LUDecomp( var A : rmatrix,
var Index : ivector,
numRows : int,
rowSwapFlag : int ) : boolean
var i, j : int
var k, iMax : int
var large, sum : real
var z, z2 : real
var scaleVect : rvector
% initialize row interchange flag
rowSwapFlag := 1
% loop to obtain the scaling element
for i := 0...numRows-1 do
large := 0
for j := 0...numRows-1 do
z2 := abs( A[i,j] )
if z2 > large then
large := z2
end if
end for % j
% no non-zero large value? then exit with an error code
if large = 0 then
return false
end if
scaleVect[i] := 1 / large
end for % i
for j := 0...numRows-1 do
for i := 0...j-1 do
sum := A[i,j]
for k := 0...i-1 do
sum := sum - A[i,k] * A[k,j]
end for
A[i,j] := sum
end for
large := 0
for i := j...numRows-1 do
sum := A[i,j]
for k := 0...j-1 do
sum := sum - A[i,k] * A[k,j]
end for
A[i,j] := sum
z := scaleVect[i] * abs( sum )
if z >= large then
large := z
iMax := i
end if
end for % i
if j ~= iMax then
for k := 0...numRows-1 do
z := A[iMax,k]
A[iMax,k] := A[j,k]
A[j,k] := z
end for % k
rowSwapFlag := -1 * rowSwapFlag
scaleVect[iMax] := scaleVect[j]
end if
Index[j] := iMax
if A[j,j] = 0 then
A[j,j] := MATVECT_EPSILON
end if
if j ~= numRows then
z := 1 / A[j,j]
for i := j + 1...numRows-1 do
A[i,j] := A[i,j] * z
end for % i
end if
end for % j
return true
end function
%
% This function uses the LU decomposition of a matrix to solve:
%
% A x = b
%
% the solution for x is returned in b.
%
procedure LUBackSubst( var A : rmatrix,
var Index : ivector,
numRows : int,
var b : rvector )
var i, j : int
var idx, k : int
var sum : real
k := -1
for i := 0...numRows-1 do
idx := Index[i]
sum := b[idx]
b[idx] := b[i]
if k > -1 then
for j := k...i-1 do
sum := sum - A[i,j] * b[j]
end for % j
elsif sum ~= 0 then
k := i
end if
b[i] := sum
end for % i
for decreasing i := numRows-1...0 do
sum := b[i]
for j := i + 1...numRows-1 do
sum := sum - A[i,j] * b[j]
end for % j
b[i] := sum / A[i,i]
end for % i
end procedure
%
% returns the absolute value of the argument
%
function abs( x : real ) : real
if x < 0.0 then
x := -x
end if
return x
end function
%
% returns a complex number as a string
%
function cstr( var c : complex ) : string
var s : string
var re, im : real
re := c.re
im := c.im
if im >= 0.0 then
s := frealstr(re,14,9 ) & " +" & frealstr(im,12,9 ) & "j"
else
im := -im
s := frealstr(re,14,9 ) & " -" & frealstr(im,12,9 ) & "j"
end if
return s
end function
